The Spatial Aggregation Language       
by Feng Zhao, Christopher Bailey-Kellogg, Ivan Ordonez

Listing One
// (1.1) Read input points.
points = read_points(infile)

// (1.2) Aggregate into a Delaunay triangulation; find the
//  minimal spanning tree.
point_trig = aggregate(points, make_mesh_delaunay())
point_mst = mst(point_trig)

// (1.3) Classify, keeping close-enough neighbors.
good_adj = function(adj) {
  nearby_ngraph = reachable_ngraph(point_mst, adj, depth)
  nearby_edges = faces(make_ngraph_geom(nearby_ngraph), 1)
  edge_lengths = map(nearby_edges, volume)
  m = average(edge_lengths)
  sigma = stddev(edge_lengths, m)
  distance(from(adj), to(adj)) < m + num_devs*sigma
}
point_classifier = classify(points,
make_classifier_transitive(point_mst,good_adj))
// (2.1) Redescribe point classes as trajectories.
points_to_traj = redescribe(classes(point_classifier),
                   make_redescribe_op_path_nline(point_mst))
trajs = high_level_objects(points_to_traj)

// (2.2) Aggregate trajectories based on point ngraph.
traj_ngraph = aggregate(trajs, make_ngraph_connected_substructure(point_mst,
                                          points_classifier, points_to_traj))
// (2.3) Classify based on similarity of tangent vectors.
same_limit = function(adj) {
  t1 = from(adj); t2 = to(adj)
  p1a = end1(t1); p1b = end2(t1)
  t2_points = localize(points_to_traj, t2)
  t2_corresponding_p1a = intersection(t2_points, neighbors(point_trig, p1a))
  t2_corresponding_p1b = intersection(t2_points, neighbors(point_trig, p1b))
  if (size(t2_corresponding_p1a) == 0 || size(t2_corresponding_p1b) == 0) {
    false
  } else {
    tan_1a = tangent(t1, p1a); tan_1b = tangent(t1, p1b)
    tan_2a = space_centroid(map(p in t2_corresponding_p1a,{ tangent(t2, p) }))
    tan_2b = space_centroid(map(p in t2_corresponding_p1b,{ tangent(t2, p) }))
    ((dot(tan_1a, tan_2a)>=angle_thresh} && dot(tan_1b,tan_2b)>=angle_thresh})
||
     (-dot(tan_1a, tan_2a)>=angle_thresh} && -dot(tan_1b,
tan_2b)>=angle_thresh}))
  }
}
traj_classifier = classify(trajs, make_classifier_transitive(traj_ngraph,
same_limit))